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Calculations of whistler ray paths in the outer ionosphere are shown for a variety of 
electron density profile models including exponential, constant, and columnar profiles. 
The Haselgrove formulation of the ray equations was used with the magneto-ionic repre- 
sentation of the wave refractive index to develop a set of differential equations for ray t racing 
suitable for inhomogeneous, anisotropic medium. The variation of paths with frequency, 
latitude, initial wave-normal angle, and other variables are examined for the purpose of 
providing a preliminary basis for comparison of the theoretical with some of the experi- 
mental results. 



1. Introduction 

Little work has previously been done to describe 
the mechanics and physics of electromagnetic wave 
propagation in an inhomogeneous anisotropic, me- 
dium. Of particular interest are certain phenomena 
taking place at frequencies below the plasma and 
gyrofrequeneies of the ionosphere. The subject 
has become important during the last few years 
because of the increasing interest in VLF wave 
phenomena such as "whistlers." 

Whistlers propagate in what is known as the 
"extraordinary" mode in magneto-ionic theory 
[Ratcliffe, 1959]. This mode of propagation is 
also described by the terms " whistler-mode" or 
"magneto-ionic duct." Details of the phenomena 
and additional references are given elsewhere [Helli- 
well and Morgan, 1953]. The basic theory of 
whistler propagation was first given by Storey 
[1953]; however, his description is valid only for a 
restricted set of conditions, many of which do 
not apply over the complete whistler path. 

The work reported here is an attempt to demon- 
strate some of the complex phenomena of whistler 
propagation. The set of computations described 
here is in no sense complete or comprehensive; it 
is meant only to demonstrate a few interesting 
details of this mode of propagation and to stimulate 
further thought and work. Very little attempt 
has been made to interpret the form and shape of 
the results shown; such interpretations must arise 
from an examination of the form of the general 
solution of the ray path equations. 

2. Discussion 
2.1. General Procedure 

Equations which are in a form useful for compu- 
tation of ray paths in the ionosphere can be obtained 
from substitution of the familiar Appleton-Hartree 
expression for the complex refractive index into 



i Based on a report prepared by Stanford Research Institute for Stanford Uni- 
versity under Prime Contract AF 18(603)-126 (Dec. 1959). 



ray equations derived by Haselgrove [1954]. The 
resulting equations in spherical coordinates r, 6, <f> are 
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/x=the real part of the complex phase 
refractive index, 
Pry pe, P0=the physical components of a vector 
of length fx, and directed normal to 
the phase fronts (also called the wave 
normal or refractive index vector), 

t=time of phase travel along the ray 
(that is, (fAt)[c = number of wave- 
lengths in the medium along the 
ray path), 

/=wave frequency, 

c=speed of light. 

The quantity ju and its derivatives are calculated 
from the Appleton-Hartree formula as follows : 
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where 

M= complex phase refractive index=/i— jk, 

X= normalized electron density=(iv6 2 )/(cd 2 ra€o), 

Y= normalized magnitude of the earth's magnetic 

field vector, Y=(p eH e )/(a)m), 

Z= normalized collision frequency=(v )/co, 

-> -» 

^=angle between p and F=cos~~ 1 [(p r Y T -r- peYe 

+ P ,YJ/(»Y)], 

N= electron density, 

e = charge on an electron, 

m=mass of an electron, 

co=27r/=2x • wave frequency, 

€ = dielectric constant of free space, 

/*o= permeability of free space, 

v = collision frequency in collisions per second, 

K= imaginary part of the complex phase refractive 
index. 

The two values of M corresponding to the plus and 
minus sign on S represent the two modes of iono- 
spheric propagation commonly called the "ordinary" 
and "extraordinary" modes. 
The derivatives of p are 

dp _dp d^_dp /p t Y cos ^—pYA 
dp~df dp~dt \ p 2 Y sin $ ) 

where t=the coordinates r, 6, </>. 

[Note: When ^-»0, d M /d^->0, b^/d Pr ->«> but d/i/d Pr ->0.] 
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where Im = imaginary part of. 



[Note: d\f//di is calculated holding p constant and so measures 
the change in the direction of the earth's magnetic field in 
space.] 

The quantities bX/di and dZ/bi must be derived 
from the space variation of X and Z which are 
assumed. An arbitrary inhomogeneous ionosphere 
can be represented in this manner. If an earth- 
centered dipole field approximation is used for the 
magnetic field, then 
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where Y e =Y at the equator at the surface of the 
earth « (8.7 • 10 5 )//. The vector components of Y 
in spherical coordinates are 
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0=colatitude and 
<p=longitude. 

The derivatives of \f/ for a dipole magnetic field are 



br 



=0 
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-> 
where 7= angle between Y and the radial (see fig. 1). 

[Note: 7= tan -1 (— }& tan 0) for the dipole field.] 
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[Note: When ^— >0, both numerator and denominator of 
dxf/jdy go to zero and bij/jby—*l.] 
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/ wave normal vector p 

magnetic field line vector Y 



It is usually desired to compute the following 
additional quantities along with the ray path: 
(1) Path length, S 

dS 



at fi z ray refractive index p cos a 

(2) Time of travel, T, also called group delay 
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2.2. Notes on the Method of Digital Computation Used 

a. Initial Conditions Needed 

Initial values for r, 6, </>, p,, p e , p are needed. To 
get p r , pe, and p^,, the wave normal direction can be 
inserted as an initial condition, p. is then calculated 
and the p's obtained. For example, let x t =(pi/p) = 
direction cosines of wave normal with respect to the 
local coordinate axes. Then, since p=p, pi=pXi. 

b. Stability Condition 

The six basic ray differential equations contain 
one more degree of freedom than there are con- 
straints explicitly expressed. That is, it only takes 
5 quantities to define a vector field (three space 
coordinates and two direction cosines). The sixth 



Figure 1. Coordinate system for two 
dimensions, r-0. 



equation is implicit in the definition of |p | =/x- The 
round-off errors involved in calculating p and in 

integrating dpi/dt will result in |p|^M after a few 
steps of integration unless an additional constraint 
is imposed. Thus, to maintain the equality, we 
correct the p/s at each point in the manner 

— ft 

Pi corrected Pi 
P 

This correction is made as early as possible in the 
derivative calculation, i.e., immediately after calcu- 
lation of M. Since the magnitude of p is not used 
in the calculation of M (only the wave normal direc- 
tion is used), no errors arising from \p\ f^/x are allowed 
to accumulate by this procedure. 

c. Interval Size 

The A/ to use depends upon the accuracy required, 
the frequency, the size of refractive index gradients 
encountered along the ray path, and the stability 
and accuracy of the numerical method of solution of 
the differential equations used. For the whistler 
ray paths, shown later in this report, the interval 
size was changed when the local truncation error 
became either too large or too small. The local 
truncation error was measured by noting the differ- 
ence between the predicted values for the next point 
and the values obtained by application of the first 
corrector (an Adam's predictor-corrector method of 
solution of the D.E. equations was used). When 
this difference became less than 5- 10~ 7 , the interval 
size was doubled using an extrapolation formula in a 
special subroutine. When the difference became 
greater than 10 ~ 4 , the interval size was halved. The 
corrector was applied only once every 5 or 10 points 
and only a predictor was used for the in-between 
points. This was done because application of a 
corrector and subsequent doubling or halving costs 
several "points" worth of time. Thus it is undesir- 
able to do this very often. The interval size used 
for these whistler calculations was approximately 
one-half an earth radius. Calculation was started 
at the top of the F 2 layer. 
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d. Other Comments 

The bulk of computation time is spent in calculat- 
ing the derivatives of M, /x, etc. Only a small frac- 
tion of the total time is actually spent in the numeri- 
cal method of solution calculations. Therefore, a 
method of solution should be chosen that minimizes 
the number of times the derivatives have to be 
calculated, and that makes maximum use of the 
derivatives already calculated. Thus the Runge- 
Kutta methods are long and cumbersome for these 
equations. A predictor-corrector method with a 
minimum number of corrector cycles applied seems 
to be about the best. Perhaps even a high-order 
predictor (5th or 6th order) might be best, since this 
would increase the interval size one could use without 
increasing the calculation time proportionately. 

The main weakness of the program used for these 
calculations is the instability of the extrapolation 
formula used for doubling. This formula has to 
jump ahead two intervals. It would be better to 
carry more points along with the calculations and 
use a higher-order extrapolation formula that only 
jumps ahead one space — or better still, carry suffi- 
cient points to use no extrapolation formula at all. 

If these equations are to be used for frequencies 
high enough that the earth's magnetic field can be 
ignored (frequencies much higher than the local 
gyrofrequency) then a separate routine should be 
written for this case, since the equations simplify to 
a much easier and faster computation. If the same 
routine is used, many of the expressions for the 
derivatives become indeterminate when Y=0. Most 
of these derivatives actually become zero, but special 
tests would have to be installed in the program to 
handle this case. 

2.3. Computations 
a. General Procedure 

The differential equations were programed for 
numerical solution on an Electrodata 205 computer. 
Only the two-dimensional case of propagation in the 
magnetic meridian was utilized in these computa- 
tions. Also, no losses were included (v=0). The 
calculation of a whistler path was begun with the 
assumed initial conditions at an altitude of 300 km. 
This altitude corresponds roughly with the maximum 
ionization level of the F 2 layer. 

An electron density profile was assumed for alti- 
tudes greater than this. 

The computations are organized into the following 
categories according to the electron density profile 
used and the parametric variation examined: 

(1) Exponential model— variation of frequency, 

(2) Exponential model — variation of initial latitude, 

(3) Exponential model — variation of initial wave 

normal direction, 

(4) Constant density model, and 

(5) Miscellaneous. 

The exponential model used in all of the first three 
categories is: 



N= 180,000 exp [-4.183119 (r- 1.0471)] 

where 

N= electron density in electrons per cubic centimeter. 
r= radial distance from center of earth in earth 
radii. 

(r= 1.0471 = 300 km above the surface of the earth 
at the beginning of the path computation). 
These particular values of electron density and scale 
height were obtained from Maeda and Kimura [1956]. 

b. Exponential Model — Variation With Frequency 

Figure 2 (a) gives the results of ray path computa- 
tions at the frequencies 5, 10, 20, 25, 30, 50, 100, 
200, 400, 1000 kc/s. Figures 2 (a) and (b) show the 
ray path for each of these frequencies. The initial 
wave normal direction for each of these paths is 
vertical (wave normal angle = 0°). In order to 
exhibit the variation of certain variables in these 
computations, figures 2(c) to2(i) have been plotted. 

Some of the interesting features of these computa- 
tions may be listed as follows: 

c. Exponential Model — Variation With Latitude 

Figure 3 gives the results of ray path computations 
initiating at north latitudes of 10, 20, 30, 35, 40, 45, 
50, and 60°. Figures 3 (a) to (i) show the ray paths 
for each of these latitudes as well as the variation of 
wave normal direction over the path. Additional 
data from these computations have been plotted in 
figures 3 (j) to (o). 

d. Exponential Model — Initial Wave Normal Variation 

Figure 4 gives the results of ray path computations 
with initial wave normal angles 10, 5, —5, and —10° 
measured positively clockwise from vertical (along 
the radial). Figures 4 (a) to (d) show the ray paths 
and wave normal direction over these paths. The 
wave frequency is 10 kc/s. Additional data from 
these computations have been plotted in figures 4 (e) 
to (g). 

The behavior of the wave normal over the path 
changes only negligibly with these changes in initial 
wave normal direction. 

e. Constant Density Model 

Figure 5 gives the results of ray path computations 
with constant electron density. The characteristics 
of a constant density model are as follows. Since 
there is no electron density gradient, the total refrac- 
tive index gradient is formed by the magnetic field. 
Also, since the refractive index gradient determines 
the change in wave normal which in turn determines 
the bending of the ray path, the bending of the ray 
path is determined only by the guiding effect of the 
magnetic field, not by the electron density. Thus 
the ray path is the same regardless of the magnitude 
of the electron density. 
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f. Miscellaneous Computations 

Comparisons with Maeda's computations — Maeda 
(1951) used a very approximate method of numerical 
integration as well as an approximation to the refrac- 
tive index gradient which is only very roughly 
correct . In order to compare his computations with 
the more exact formulation of Haselgrove, one of 
Maeda 's published ray paths was compared with a 
similar computation using the present program in 
figure 6. It can be seen that the difference is 
appreciable. The lack of comparison is not due to the 
quasi-longitudinal approximation of the refractive 
index function which was used by Maeda since this 
approximation was used also in Haselgrove's formula- 
tion for this particular computation. The difference 
between the Q. L. approximation and the full expres- 
sion is too small to be shown. The differences 
between the two computed paths may be due to the 
relatively crude method of integration used by Maeda 
to numerically solve the differential equations. 

Field-alined column model— In a recent theoretical 
paper it lias been shown that if the ionization is 
field-alined it can act like a waveguide in trapping 
whistler energy [Smith, Helliwell, and Yabroff, I960]. 
This would explain some of the very distinct and 
clear whistlers which have been recorded. In order 
to answer certain questions regarding the effects of 
field -alinement, the following models have been 
computed: 



iV=180,000/l+Cexp P^T"]} 

exp [-4.183119(r-1.0471)] 
where 

b = r/sin 2 d (the equation of a field line is b== 

constant) , 
b = 2.094 (this is a field line corresponding to 300 
km altitude at a latitude of 45° and is the 
center of the column), 
(7= the modulation factor, i.e., the relative 
increase of maxim inn ionization at the 
center of the column over that of the sur- 
rounding background level of ionization, 
and 
D=ihe standard deviation of the column. 
(Tl ie column is of gaussian electron density 
distribution with the standard deviation 
measured in terms of the local distance 
between the field lines.) This makes more 
intuitive physical sense than making the 
thickness of the column constant at all 
points on the field line, since the individual 
electrons will tend to diffuse easily along 
the lines of force but not across them. 
For these computations D=0.029 which corres- 
ponds to a standard deviation of about 35 km at an 
altitude of 300 km, and about 215 km at. the top of 
the path. 




Figure 2(a). Ray paths at various 
frequencies. 
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INITIAL WAVE NORMAL ANGLE * 0' 



Figure 2(b). Ray paths at various 
frequencies. 
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Figure 2(c). Maximum altitude versus frequency. 

The maximum height of the path remains essentially constant up to the frequency of minimum time 'delay [the nose 
frequency, which is seen to be at about 20 kc/s in fig. 2(d)] at which point the maximum path height begins decreasing. 
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INITIAL LAT ITUDE ■ 45 N 
INITIAL WAVE NORMAL ANGLE 



FREQUENCY, kc/s 



Figure 2(d). Total path delay versus frequency. 



The time delay decreases in the expected "Eckersley" fashion until the nose region is reached. This region corresponds Jto 
the frequency reaching approximately 45 percent of the gyrofrequency at the top of the path. These calculations confirm this 
behavior. The time delay at higher frequencies Increases until the elTect of the shortening of the path length at high frequencies 
overbalances the decrease in group velocity to cause the time delay to begin to decrease again above 200 kc/s. 
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Figure 2(e). Final latitude versus frequency. 

The final latitude remains fairly constant at frequencies up to the nose frequency. This suggests some latitude focusing in 
the "Eckersley" region. Above the nose frequency, the final latitude starts shifting toward the pole until the effect of the 
shortening path length at higher frequencies overcomes this tendency and results in a shift back toward the equator. 
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Figure 2(f). Wave normal — field line angle along the path. 

At frequencies below the nose frequency, the wave normal swings steadily 
below the field line as the path progresses. As the nose frequency is approached, 
the wave normal tends to oscillate about a position approximately ten degrees 
above the field line. At higher frequencies, the wave normal swings outward 
after an initial dip toward the field line. 
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Figure 2(g). Wave normal — field line angle along the path. 
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Figure 2(h) . Wave normal — field line angle 
along the path. 



In the first two computations shown in figures 7 
and 8 the initial wave normal direction is vertical. 
When the modulation is 0.2, good guiding occurs. 
When the modulation is reduced to 0.1, the ray is not 
guided because the initial wave normal direction is 
far from the field-line direction and the ray gets 
pulled out of the column before the wave normal gets 
properly lined up with the field line. This is in 
agreement with theoretical predictions [Smith, 
Helliwell, and Yabroff, I960]. In the third case 
shown in figure 9, the wave normal was set along the 



field line. In that case, the guiding is even better 
than for the higher modulation of case 1. The final 
case shown in figure 10 shows the modulation reduced 
to zero with the initial wave normal direction on the 
field line. The guiding of this ray is still good, but 
less so than for the previous case of 0.1 modulation. 
These calculations support the theoretical pre- 
dictions [Smith, Helliwell, and Yabroff, 1960] that 
wave components whose wave normals lie within a 
certain angle with respect to the direction of the 
magnetic field will be trapped. 
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Figure 2(i). Wave normal — field line angle 
c long the path. 



Figure 3(a). Ray path initiating at 10 °N 
latitude. 
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f = 5kc/s 

INITIAL WAVE NORMAL ANGLE ■ 0' 

INITIAL LATITUDE = 20° N 



Figure 3(b). Ray path initiating at 20 °N 
latitude. 



Figure 3(c). Ray path initiating at 30 °N 
latitude. 
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f=5kc/s 
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Figure 3(d). Ray path initiating at 35 °N 
latitude. 
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Figure 3(e). Ray path initiating at /+0 °N 
latitude. 
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f = 5kc/s 
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Figure 3(f). Ray path initiating at 45 °N 
latitude. 



Figure 3(g). 



Ray path initiating at 50 °N 
latitude* 
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Figure 3(h) . Ray path initiating at 55° N 
latitude. 
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Figure 3(i). Ray path initiating at 60 °N so'< 
latitude. 
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f« 5kc/s 

INITIAL WAVE NORMAL ANGLE = 0° 



Figure 3(j). 



INITIAL NORTH LATITUDE, deg 

Maximum altitude versus initial latitude. 



The maximum height of the path is seen to increase very rapidly with increas- 
ing initial latitude until an altitude is reached at which the plasma frequency 
becomes equal to the wave frequency. This represents the usual reflection level. 




INITIAL NORTH LATITUDE , deg 

Figure 3(k). Total path delay versus initial latitude. 

The total time delay is seen to be fairly constant over a wide range of initial 
latitude. This may be a bit surprising since it was previously observed that the 
path length changes very rapidly with initial latitude. Clearly the average 
group velocity must be much higher for the longer paths which are at a higher 
latitude. 
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Figure 3(1). Final latitude versus initial latitude. 

This plot of final latitude against initial value shows some focusing of rays in 
the region of 45° and 33° latitude. 
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Figure 3(m). Wave normal — field line angle along the path. 

At initial north latitudes greater than 35°, the wave normal swings from the 
outside to the inside of the field line along the path. For initial north latitudes 
less than 35°, the wave normal swings farther away from the field line. 
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Figure 3(n). Wave normal — field line angle along the palh. 
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Figure 3(o). Wave normal — field line angle along Ike path. 



Figure 4(a). Ray path with initial wave 
normal angle= 10° . 
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Figure 4(b). Ray path with initial wave 
normal angle=5° 



Figure 4(c). Ray path with initial wave 
normal angle— —5°. 
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Figure 4(d). Ray path with initial wave 
normal angle = 10°. 
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Figure 4(e). Maximum altitude versus initial wave normal 
angle. 

Maximum path height is seen to decrease with increasing wave normal angle 
until the wave reaches between 5° and 10°. 
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Figure 4(f). Total path delay versus initial wave normal 

angle. 

The total time delay is also seen to remain roughly constant with changes in 
initial wave normal angle. This is again consistent with the idea of average 
group velocity being roughly proportional to path length. 
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Figure 4(g). Final latitude versus initial wave normal angle. 

The final latitude is seen to remain constant with initial wave normal variation 
even though maximum path height changes. 



Figure 5. Constant density model. 
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Figure 6. Comparison with Maeda's 
computation. 



Figure 7. A field-alined column model. 
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Figure 8. A field-alined column model. 



Figure 9, A field-alined column model. 
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Figure 10. Ray path with initial wave 
normal along field line. 
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